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i^\BSTRACT: Strain in granular materials in quasistatic conditions under varying stress originate in (I) contact 
Qleformation and (II) rearrangements of the contact network. Depending on sample history and applied load, ei- 
ther mechanism might dominate. One may thus define rheological regimes I and II accordingly. Their properties 
e presented and illustrated here with discrete numerical simulation results on sphere packings. Understanding 
he microscopic physical origin of strain enables one to clarify such issues as the existence of macroscopic 
L^-elasticity, the approach to stress-strain relations in the large system limit and the sensitivity to noise. 



«tt INTRODUCTION 

^(acroscopic strain in solidlike granular materials has 
two obvious physical origins: first, grains deform near 
their contacts, where stresses concentrate (so that one 
yodels intergranular interaction with a point force); 
TZthen, grain packs rearrange as contact networks, be- 
tween two different equilibrium states break, and then 
Ojepair in a different stable configuration. We refer here 
1 — respectively to the two different kinds of strains as 
type I and II. The purpose of the present communi- 
cation is to delineate the regimes, denoted as I and 
*Ol accordingly, within which one mechanism or the 
^yther dominates, in a simple model material (an as- 
sembly of spheres), from discrete numerical simula- 
tions. Macroscopic mechanical properties are shown 
Qd differ, as well as microscopic variables. 

| The very small strain elastic response of granu- 
lar materials belongs to regime I: what is measured 
. £jhen is the macroscopic stiffness of a spring net- 
^vork, each intergranular contact behaving like an 
plastic element. Such a spring network model is usu- 
ally ado pted on studying vibration modes and elasti c 
moduli (ISomfai et al. 20041 |Agnolin & Roux 2005[ ). 
However, elastic-frictional contact networks also 
comprise plastic elements (sliders), and deform ir- 
reversibly under quasistatically applied stress incre- 
ments. As long as they still support the applied load, 
strain amplitudes scale as the inverse of the stiffness 
constants of the springs. Such a scaling will be used 
here as a signature of regime I, which extends, be- 
yond the quasi-elastic domain, throughout the stress 
or strain interval corresponding to the elastoplastic re- 
sponse of a given contact network. If grains are mod- 
eled as perfectly rigid, strains in regime I reduce to 
zero. 

Regime II will in general correspond to larger 



strains, for which contact networks keep rearrang- 
ing. Strain amplitudes are then related to the distances 
(gaps) between neighbouring grains that do not touch. 
Contact stiffnesses are then expected to have little 
influence on macroscopic deformations. Such situa- 
tions are sometimes studied by simulation methods 
that deal with rigid grains, such as Contact Dynamics 
(Radjai & Roux 2004). As the network continuously 
fails and repairs, larger dynamical effects and larger 
spatial fluctuations of strain are expected, since fail- 
ing materials usually exhibit larger heterogeneities. 

The simulations reported below explore the condi- 
tions of occurrence of regimes I and II, and give a 
quantitative meaning to the statements made in this 
introduction. After basic features of numerical com- 
putations are introduced in Sec.[2l results on the con- 
stitutive law and regimes I and II are given in Sec. [3l 
Sec.[4]is a brief conclusion. 

2 NUMERICAL MODEL 

Triaxial compression tests of monosized assemblies 
of N (N = 4000 for most results here), spheres of 
diameter a were simulated by molecular dynamics 
(MD, or DEM). In those computer experiments, one 
starts with an isotropically assembled initial state with 
pressure P, and then, keeping the axes of coordinates 
as principal stress directions, increases slowly the 
largest principal stress, a±, while the others are held 
fixed, equal to P. One denotes as q the stress deviator, 
q = o"i — p. Like in most numerical studies, we chose 
here to impose a constant strain rate i\ and to mea- 
sure o\ as a function of e\, termed "axial strain" and 
subsequantly denoted as e a . Soil mechanics conven- 
tions are adopted: compressive stresses and shrinking 
strains are positive. We focus on the quasistatic me- 
chanical behaviour expressed by dependences q(e a ), 
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e v (e a ) as e a increases, in dense systems, before the de- 
viator peak is reached. t v is the volumetric strain (rel- 
ative volume decrease). Dimensional analysis leads to 
the definition of the inertia parameter I = e a ^Jm/aP 
as a measure of the departure from equilibrium, the 
quasistatic limit being I — > 0. 

Motivated by possible comparisons to laboratory 
experiments with glass bead packings, simulations 
are carried out with Hertz-Mindlin contacts, with 
the elastic properties of glass (Young modulus E = 
70 GPa, Poisson coefficient v = 0.3), and a friction 
coefficient fi = 0.3 - additional details and references 



are provided in ( |Agnolin & Roux 2005 [ ). Normal vis- 
cous forces are also implemented: the damping pa- 
rameter in any contact is chosen as a fixed fraction 
C of its critical value, defined for the contacting pair 
with its instantaneous (i.e., dependent on current nor- 
mal force or h) stiffness constant dF^/dh. A suit- 
able dimensionless parameter characterizing the im- 
portance of elastic deflections h in contacts is k = 
(E/P) 2 / 3 (such that h/a ~ 1/k). 

3 SIMULATION RESULTS 

3.1 Sample preparation 

In order to obtain dense samples, we first simulated 
sets of sphere packings prepared by isotropic com- 
pression of frictionless granular gases. This results in 
configurations hereafter denoted as A. A-type config- 
urations have a high coordination number (approach- 
ing 6 at low pressure if inactive grains are discarded). 
They therefore present a large force indeterminacy. 
We observed (fig. Q]) that the raise of deviator q with 
axial strain in such samples is much faster than in 
usual experimental results, for which e a is usually of 
order 1% to 5% at the deviator peak. Likewise, the 
onset of dilatancy after the initial contractant strain 
interval is unusually fast in A samples. This moti- 
vated the use of a different preparation procedure, 
which, although idealized, aims to imitate the effects 
of vibrations in the assembling of a dense, dry pack- 
ing of beads. In this method (called C in the se- 
quel), A samples are first dilated (multiplying coor- 
dinates by 1.005), then mixed, as by thermal agita- 
tion, until each grain has had 50 collisions on aver- 
age, and finally compressed in the presence of fric- 
tion to a relativy low pressure, P = lOkPa. Higher 
P values are obtained on further compressing. Fig. \T\ 
compares the behaviour of initial states A and C, in 
triaxial compression with P = 100 kPa (k ~ 6000). 



Agn olin & Roux (2005 [ ) report in these proceedings 
on the large difference in coordination number be- 
tween A and C states, where it is much smaller (~ 
4.7), while densities are very close. Usual experimen- 
tal curves, which do not exhibit q maxima or dila- 
tancy before e a ~ 0.01, are better modelled with C 
samples. Those experiments are made with, e.g., dry 
grains assembled in the laboratory. One cannot ex- 
clude, however, that samples left to age and anneal 
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Figure 1. q(e a ) (left scale) and e v (e a ) (right scale) curves for 
A and C states under P = 100 kPa. Averages over 5 samples of 
4000 spherical grains. 

for a long time gradually evolve towards better coor- 
dinated configurations resembling A ones. One may 
also assemble the grains in the presence of a lubri- 
cant, thereby strongly reducing friction in the initial 
stage (Agnolin et al. 2005). A-type samples can thus 
be viewed as ideal models for preparation procedures 
suppressing friction, while C ones are more appropri- 
ate models for laboratory specimens made by pour- 
ing, vibrating or tapping. A similar conclusion was 



reached by |Agnolin et al. (20 05 ) in a study of sound 
propagation velocities. 

3.2 Reproducibility, quasistatic limit 
Stress-strain curves as displayed on fig. Q] should 
express a macroscopic, quasistatic constitutive law. 
Sample to sample fluctuations should regress in the 
large system limit, and the results should be indepen- 
dent on dynamical parameters such as inertia, viscous 
dissipation, and strain rate, summarized in dimension- 
less parameters £ and /. Fig. [2] is indicative of sam- 
ple to sample fluctuations with 4000 beads. One may 
notice the very good reproducibility of the curve be- 
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Figure 2. Detail of small strain part of q(e a ) curves for 5 different 
samples of each type, A (top curves) and C (bottom ones) with 
N = 4000 beads. 

tween A samples in the initial fastly growing part. We 
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checked that differences between samples decreased 
for increasing N. As to the influence of dynamical pa- 
rameters, fig. [3] shows that the quasistatic limit is cor- 
rectly approached for / < 1CT 3 , a quite satisfactory 
result, given that usual laboratory tests with e a ~ 1CT 5 
correspond to I < 1CT 8 . In previous 2D simulations 



1.5 



0.5 



1 1 1 1 | 1 1 1 1 | 1 1 1 
Effects of dampin 

and strain rate 

yit=0.3 ; C (1 




P = cr 2 =CT 3 = lOOkPa 



0.03 



0.02 



- 0.01 



0.01 0.02 0.03 0.04 0.05 
e a 

Figure 3. Effect of dynamical parameters: q(e a ) and e v (e a ) 
curves for the different values of ( and / indicated coincide, 
showing the innocuousness of dynamical parameter choice. 



with disks (IRoux & Combe 2002|) . sample to sample 
fluctuations were shown to regress as N~ l l 2 . 

3.3 Influence of contact stiffness 

The small strain (say e a < 5.1CT 4 ) interval for A sam- 
ples, with its fast q increase, is in fact in regime I. This 
is readily checked on changing the confining pressure. 
Fig. H] shows the curves for triaxial compressions at 
different P values (separated by a factor vTO) from 
10 kPa to 1 MPa, with a rescaling of the strains by 
the stiffness parameter k, in one A sample. Their co- 
incidence for q/P < 0.8P evidences a wide deviator 
range in regime I. For larger strains, curves separate 
on this scale, and tend to collapse together if q/P, e v 
are simply plotted versus e a . The strain dependence 
on stress ratio is independent from contact stiffness. 
This different sensitivity to pressure is characteristic 
of regime II. Fig [5] shows that it applies to C sam- 
ples almost throughout the investigated range, down 
to small deviators (a behaviour closer to usual experi- 
mental results than type A configurations). At the ori- 
gin (close to the initial isotropic state, see inset on 
fig. [5]), the tangent to the curve is given by the elastic 
(Young) modulus of the granular material, E m , which 
scales as k, but curves quickly depart from this be- 
haviour (below q = 0.1P). 

3 .4 Load reversal 

If (fig. [6]) one reverses the direction of load incre- 
ments, the stress-strain curves exhibit notable in- 
tervals within which the deviator stress decreases 



0.5 



0.5 




-0.5 



0.001 0.002 0.003 0.004 
e a (P /P) 2/3 

Figure 4. q(e a )/P and e v (e a ) curves for one A sample and 
different P values. Strains on scale (P/Pq) 2 / 3 oc Pq = 
100 kPa. 
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Figure 5. q(e a )/p for same p values as on fig.|U in one C sam- 
ple. Inset: detail of same curves, blown-up e scale, straight lines 
corresponding to Young moduli in isotropic state. 

very fast, which results in large irreversible (plas- 
tic) strains. It can be checked that the initial slope of 
those descending curves are equal to the Young mod- 
ulus E m of the material, and that subsequent strains 
scale as 1/E m , like the initial q increase in A sam- 
ples. Therefore, some significant deviator stress inter- 
vals (of order 0.2P or larger) are found in regime I on 
reversing deviator stress or axial strain variations. 

Fig. U\ shows that the small strain response of A 
samples, within regime I, close to the initial state, is 
already irreversible. Type I strains are not elastic. 

3.5 Calculations with a fixed contact list 

Within regime I, the mechanical properties of the ma- 
terial can be successfully predicted on studying the 
response of one given set of contacts. Those might 
slide or open, but the very few new contacts that are 
created can be neglected. To check this in simulations, 
one may restrict at each time step the search for inter- 
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Figure 6. Top plot: effects of load reversals at different points on 
curves (C sample). Initial slopes of unloading curves correspond 
to elastic moduli. Bottom: evolution with e a of some elastic mod- 
uli, probing induced anisotropy. 
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Figure 7. Very small strain part of q(e a ) curve in one A sample, 
showing beginning of unloading curves (arrows). Curve marked 
NCC was obtained on calculating the evolution of the same sam- 
ple without any contact creation. 

acting grains to the list of initially contacting pairs. 
Fig. |7] compares such a procedure to the complete 
calculation. The curve marked "NCC" for no contact 
creation is indistinguishable from the other one for 
q > 0.8. In two dimensions, |Roux & Combe (2002[ ) 
could implement a purely static method (elastoplastic 
computation on a given contact network), apt to cal- 
culate the quasistatic evolution of the sample under 
varying applied stresses throughout the initial regime 
I stage of 2D assemblies of disks analogous to A sam- 
ples. The limit between regimes I and II was studied 
with some accuracy (|Combe 2002b . and shown to ap- 
proach a finite value in the rigid limit (k — > +oo), and 
in the limit of large systems. This value does not ap- 
pear to depend on details of contact elast icity, suc h as 
tangential to normal stiffness ratio (|Combe 2002|) . 

3.6 Microscopic aspects 

The existence of wide stress intervals within regime 
I is associated with strongly hyperstatic contact net- 
works (large force indeterminacy). Initially, A sam- 



( Agnolin & Roux 2005 ), and friction is not mobilized 
(zero tangential forces). Consequently, the set of con- 
tact forces that resolve the load and satisfy Coulomb 
inequalities is large, and the initial forces are far from 
its boundaries. At coordination 6 this set spans an 
affine space of dimension 3iV* if N* is the number 
of force-carrying particles. In regime II, regarding the 
Coulomb condition in sliding contacts as a constraint 
on force values in the count of force indeterminacy, 
this dimension decreases to a fraction of order 10% 
of the number of degrees of freedom. Upon revers- 
ing the load variation, sliding contacts tend to dis- 
appear, leading to a larger force indeterminacy and 
a notable type I interval. The small variation of co- 
ordination number in the pressure range of Fig. [5] 
( |Agnolin et al. 2005] ) witnesses the smallness of ge- 
ometrical changes, hence the collapse of curves with 
type II strains. 

Larger strain heterogeneities and sensitivity to per- 
turb ations are oth er characteristic features of regime 
II (IRoux & Combe 2003D . 

4 CONCLUSION 

Numerical studies thus reveal that the two regimes, in 
which the origins of strain differ, exhibit contrasting 
properties. On attempting to predict a macroscopic 
mechanical response from packing geometry and con- 
tact laws, the information about which kind of strain 
should dominate is crucial. Regime I corresponds in 
usual testing conditions to highly coordinated systems 
(with many contacts), or to changes in the direction 
of load increments (hence a loss in friction mobiliza- 
tion). Investigating the nature of strains might open 
interesting perspectives to study the effects of cyclic 
loadings or random perturbations. 
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